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Multiresolution  Representation  and  Numerical  Algorithms: 

A  Brief  Review 

Ami  Harten1 

School  of  Mathematical  Sciences 
Tel- Aviv  University 
Tel- Aviv,  69978  ISRAEL 
and 

Department  of  Mathematics 
UCLA 

ABSTRACT 

In  this  paper  we  review  recent  developments  in  techniques  to  represent  data  in  terms  of  its 
local  scale  components.  These  techniques  enable  us  to  obtain  data  compression  by  eliminating 
scale-coefficients  which  are  sufficiently  small.  This  capability  for  data  compression  can  be 
used  to  reduce  the  cost  of  many  numerical  solution  algorithms  by  either  applying  it  to  the 
numerical  solution  operator  in  order  to  get  an  approximate  sparse  representation,  or  by 
applying  it  to  the  numerical  solution  itself  in  order  to  reduce  the  number  of  quantities  that 
need  to  be  computed. 
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Contract  Nos.  NAS1-18605  and  NAS1-19480  while  the  author  was  in  residence  at  the  Institute  for  Computer 
Applications  in  Science  and  Engineering  (ICASE),  M/S  132C,  NASA  Langley  Research  Center,  Hampton, 
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1.  Introduction 


Fourier  analysis,  which  provides  a  way  to  represent  square-integrable  functions  in  terms 
of  their  sinusoidal  scale-components,  has  contributed  greatly  to  all  fields  of  science.  The 
main  drawback  of  Fourier  analysis  is  in  its  globality  —  a  single  irregularity  in  the  function 
dominates  the  behavior  of  the  scale-coefficients  and  prevents  us  from  getting  immediate 
information  about  the  behavior  of  the  function  elsewhere. 

The  recent  development  of  the  theory  of  wavelets  (see  [Me]  and  [Ma])  was  a  great  step 
towards  local  scale  decomposition,  and  has  already  had  great  impact  on  several  fields  of 
science.  In  numerical  analysis  representation  by  compactly  supported  wavelets  (see  [Da]  and 
[CDF])  is  used  to  reduce  the  cost  of  many  numerical  solution  algorithms  by  either  applying  it 
to  the  numerical  solution  operator  to  obtain  an  approximate  sparse  form  (see  [BCR]),  or  by 
applying  it  to  the  numerical  solution  itself  to  obtain  an  approximate  reduced  representation 
in  order  to  solve  for  less  quantities  (see  e.g.  [LT],  [MR]  and  [BMP]).  The  main  drawback  of 
the  theory  of  wavelets  is  that  it  attempts  to  decompose  any  square  integrable  function  into 
scale-components  which  are  translates  and  dilates  of  a  single  function.  Consequently  there 
are  conceptual  difficulties  in  extending  wavelets  to  bounded  domains  and  general  geometries. 
Furthermore,  the  uniformity  of  the  underlying  wavelet  approximation  makes  it  impossible  to 
obtain  an  adaptive  (data-dependent)  multiresolution  representation  which  fits  the  approx¬ 
imation  to  the  local  nature  of  the  data.  The  only  adaptivity  which  is  possible  within  the 
theory  of  wavelets  is  through  redundant  “dictionaries.” 

In  a  series  of  works  [Hl-3]  we  have  studied  the  question  of  how  to  represent  discrete  data 
which  originates  from  unstructured  grids  in  bounded  domains  in  terms  of  scale  decompo¬ 
sition.  Combining  ideas  from  multigrid  methods,  numerical  solution  of  conservation  laws, 
hierarchial  bases  of  finite  element  spaces,  subdivision  schemes  of  Computer-Aided  Design 
and  of  course  -  the  theory  of  wavelets,  we  came  up  with  the  more  general  concept  of  “nested 
sequence  of  discretization.”  Given  discrete  data  which  can  be  associated  with  a  nested 
sequence  of  discretization  we  show  that  it  has  a  multiresolution  representation,  i.e.,  a  one- 
to-one  correspondence  between  the  given  data  and  its  scale-decomposition.  This  framework 
is  a  generalization  of  the  theory  of  wavelets  in  the  sense  that  under  conditions  of  uniformity 
its  natural  result  is  wavelets. 

In  this  paper  we  review  the  work  in  [HY],  [AC],  [ACD],  [H4-5]  where  the  previously  men¬ 
tioned  works  on  numerical  solution  algorithms  with  representation  by  wavelets  are  extended 
to  the  more  general  framework  of  nested  discretization. 
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2.  Nested  Discretization 


In  this  section  we  describe  the  class  of  discrete  data  for  which  we  can  obtain  representation 
in  terms  of  a  scale-decomposition.  We  start  with  two  examples. 

Example  1. 

Let  us  consider  continuous  functions  /  in  the  interval  [0,1] 

fef  =  c°[  o,i], 

and  let  {Xk}%=0  be  the  following  nested  sequence  of  uniform  grids 

X*  =  {xf}i0,  xk  =  ihk,  hk  =  2~kh0,  Jk  =  2* Jo  (2.1a) 

for  some  integer  J0  with  h0  =  1/Jo ■  Here  k  =  L  is  the  finest  grid  and  k  =  0  is  the  coarsest. 
Observe  that  Xk  is  obtained  by  dyadic  refinement  of  Xk_1,  i.e. 

x2i  =  s?_1,  *  =  0,...,  Jjt-i  (2.16) 

X2i-1  =  2 (^l-l  xi  *  =  1>  •  •  •  »  Jk—  \  •  (2.1c) 

Let  us  consider  now  the  discrete  values  vk  =  which  are  obtained  from  point  value 

discretization  of  f  e  T  in  the  fc-th  grid 

v*  =■■  V>kf)i  =  /(*?),  i  =  0, . . . ,  Jk,  0  <  k  <  L.  (2.2 a) 

It  follows  from  (2.1b)  that  vk~l  can  be  obtained  from  vk  by  the  decimation 


v-  1  =  i  =  0,...,  Jfc_i. 


(2.26) 


Let  Ik  j)  denote  any  continuous  function  in  [0, 1]  which  interpolates  vk~x  at  the 

grid  points  of  XA'-1,  i.e. 

/fc-1(*;t;fe-1)€^=Co[0,l]  (2.3a) 

Ik-1(xk-1-,vk-i)  =  vk-1  for  i  =  0,...,Jk-1  (2.36) 

e.g.  we  can  take  Ik~1(x;vk~1)  to  be  the  piecewise-linear  interpolation 

Ik~1(x;  vk_1)  =  tfc}  +  -  vt?),  for  x^  <  *  <  (2.4) 


Given  vk  1  we  can  approximate  vk  by 


‘(^f;^1),  i  =  0 
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Let  us  denote 


ek  =  vk  -  Ik~x(xk;  vk-x),  i  =  0,...,Jk 


and  refer  to  ek  =  {ek}i=o  as  the  prediction  error.  We  observe  that 

ek{  =  0  for  i  =  0, . . . ,  Jk-i 

and  denote  the  interpolation  error  at  the  odd  grid  points  of  Xk  by  dk  =  {dk  Jyij1 
dk  =:  ek2j_t  =  4_i  -  Ik-\xk3_ i;  vk~%  j  =  1, . .  • ,  Jk- 1. 


(2.5a) 


(2.56) 


(2.5c) 


Clearly  there  is  a  one-to-one  correspondence  between  and  {dk,vk  *}:  Knowing  vk  we  get 
vk~x  and  dk  by  (2.2b)  and  (2.5c),  respectively.  From  {dk,vk~1}  we  recover  vk  by 

f  *4  =  »f-1»  i  =  0,..., Jfc-i 

(2.6) 

l  *4-i  =  lk~l  (x2i-i ;  v*-1 )  +  4 ,  *  =  l , . . . ,  Jfc-i . 

Since  vk~x  can  be  likewise  represented  by  {dk~1 ,  vk~2}  we  get  that  there  is  a  one-to-one 
correspondence  between  the  values  of  the  finest  level  vL  and  the  sequence  of  (Jl  +  1)  elements 
{dL, . . . ,  d1,  v0}  which  we  denote  by  vm'- 

vL  <^4  {dL,...,dx,v0}  =:  vm-  (2-7) 

We  refer  to  vm  as  the  multiresolution  representation  of  vL.  It  follows  from  (2.2b)  and  (2.5c) 
that  the  direct  multiresolution  representation  (MR)  transform  vL  i— ►  vm  can  be  expressed 
by  the  following  algorithm 


DO  k  =  L, . . . ,  1 


v-  1  =  vh  i  =  0,...,  Jjt-i 

d)  =  ujj-i  -  7fc_1(®2j-i;  uA:_1)>  j  ~  •  •  •  >  •fr-v 


(2.8a) 


We  denote  this  transform  by  M,  i.e. 


vm  =  M  •  v1 


(2.86) 


Similarly  it  follows  from  (2.6)  that  the  inverse  MR  transform  vM  *-*  vL  can  be  expressed  by 


DO  k  =  l,...,L 


vki  =  vk  \  i  =  0,...,Jk- 


(2.9a) 


vkj_x-Ik  l{xkj-i,vk  1)  +  dkj ,  j  =  1,..., Jk-u 
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and  we  denote 


vL  =  M  1  •  vM.  (2.96) 

We  refer  to  dk  as  the  scale-coefficients  of  the  &-th  level  of  resolution.  For  the  piecewise- 
linear  interpolation  (2.4)  we  get  from  (2.5c)  that 

dj  =  4-i  -  +  wi_1);  (2-lOa) 

in  terms  of  /  G  T  for  which  vk  —  Vkf  this  can  be  expressed  by 

<$(/)  =  -\\f(4,)  -  2/(*5i-i)  +  .,)]•  (2.104) 

Hence  if  f(x)  is  twice  differentiable  in  we  get  that 

dj{f)  =  ~^(hk)2f"(0  for  some  (  €  (4-i»4_1)  (2.10c) 

and  consequently  the  scale  coefficients  in  a  region  of  smoothness  of  f(x)  tend  to  zero  as 
0((hk)2)  for  k  — >  oo.  However  at  a  jump  discontinuity  dk(f )  is  proportional  for  the  size  of 
the  jump  and  thus  remains  0(1)  as  k  — >  oo. 

We  can  obtain  data  compression  by  setting  to  zero  all  scale  coefficients  which  fall  below 
a  prescribed  tolerance.  Let  us  denote 


dk  =  < 

j'  0  if  \dk\  <  ek 

1 

(2.11 a) 

{  dk  if  \dk\  >  £k, 

and 

O 

£ 

\ 

II 

<s> 

(2.116) 

Based  on  the  analysis  in  [H3]  we  get  the  following  bound  on 

the  compression  error  in  the 

case  of  piecewise-linear  interpolation: 

max 

0<t<JL 

14-41  < 

k=i 

(2.12a) 

Given  any  e  >  0  we  can  take 

nk—L—l 

fit  =  2  £ 

(2.126) 

and  thus  ensure  by  (2.12a)  that 

max  {vf'  —  uf'l  <  e, 

0<i<JL  *  1 

(2.12c) 
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and  let 

C‘ =  {«?}£„  c{  =  (*L.*i).  (2.13a) 

where  {x*}  are  the  gridpoints  of  Xk  in  (2.1);  observe  that 

=  4i-i  U  4-  (2-136) 

Let  =  {vk}i± i  be  discrete  values  which  are  obtained  by  taking  the  average  over  the  cells 
in  Ck  of  some  function  /  €  T  =  L1  [0, 1] 

vk  =:  (!>*/),  =  i  A  fdx,  \ck\  =  hk,  i  =  l,...,Jk,  l<k<L.  (2.14a) 

lci  I  Jci 

It  follows  from  the  additivity  of  the  integral  and  (2.13b)  that  can  be  obtained  from  vk 
by  the  decimation 

V?'1  =  +  4),  *  =  1.  •  ■  ■ ,  Jk-1-  (2.146) 

Let  (TZkvk)(x)  denote  any  function  in  L1^,  1]  which  satisfies 

(vk •  nkvk)i  =  Jck(nkvk)(x)dx  =  vk,  i  =  i,...,Jk  (2.15) 

e.g.  we  can  take  7Zkvk  to  be  the  piecewise-constant  function 

( Kkvk)(x )  =  vk  for  x  €  ck.  (2.16) 

In  the  context  of  ENO  schemes  for  the  numerical  solution  of  conservation  laws  we  refer  to 
Kk  as  reconstruction  from  cell- averages  and  to  (2.15)  as  “conservation”  (see  [HEOC]). 
Given  vk~x  we  can  approximate  vk  by 

vk  «  (Vk  •  7efc-ivfc-1)i  =  i  [.{Kk-1  Vk~l)(x)dx, 

\ci  I  Jci 

i.e.  by  first  approximating  /  from  vk_1  by  'R.k-\vk~*  and  then  taking  cell-averages  of  this 
approximation  over  the  finer  level  k.  It  follows  from  the  conservation  property  (2.15)  that 
the  prediction  error  ek, 

ei  =  vi  —  (T^k  •  ‘ft-fc-iv*-1)*  ,  i  =  l,...,J*  (2.17a) 
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satisfies  the  relation 

e2i-i  +  4  =  0  f°r  *  =  1, 

•  •  • ,  «4_i. 

(2.176) 

Therefore,  if  we  store 

4  =  4-i  f°r  i  =  l,-- 

•  5 

(2.18a) 

we  can  recover  the  prediction  error  ek  by 

pA:  __  JA: 

e2t-l  —  at 

◄ 

(2.186) 

A  —  —dk 

e2i  ~  ai 

and  thus  get  a  one-to-one  correspondence  between  vk  and  {d*,*;*-1}.  As  in  the  previous 
example  this  leads  to  the  multiresolution  representation 


vL  e— ->  {dL,...,d\v0}  =:  vM 


where  the  direct  MR  transform  vm  —  M  •  vl  can  be  expressed  by  the  algorithm 

DO  k  =  L, . . . ,  1 

Vi  1  “  2  (4-1  +  V2i)i  *  =  1,  •  •  •  ,  Jk- 1 
dj  =  v2j- 1  -  fc*J_t('%k-ivk-1)(x)dx,  j  =  1,  .  .  .  ,Jk-i 

and  the  inverse  MR  transform  vL  =  M~l  •  vM  is  given  by 


DO  Jc  =  l,...,L 
DOi  =  l,...,Jk_1 


4-i  =  1)(x)dx-\-d^ 


k  _  Q„,A:-1 


=  2vt- 


—  v 


k 

2t— 1  • 


(2.19) 


(2.20) 


Observe  that  the  last  statement  of  (2.20)  is  obtained  from  (2.14b). 

In  [HEOC]  we  showed  that  any  interpolation  method  (2.3)  gives  rise  to  a  reconstruction 
from  cell-averages  (2.15)  by  the  following  “reconstruction  via  primitive  function”  technique: 
Given  cell-averages  vk  =  Vkf  in  (2.14a)  we  calculate 


ft  =  nd),  cm  =  I’rndy 

Jo 

by 

=  0,  Fk  =  hkJ2  vk,  1  <  i  <  Jk,  (2.21  a) 

j- 1 
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and  define 


(Rkvk)(x)  =  £lk(x-,Fk),  (2-21  b) 

where  Ik(x]  Fk)  is  any  interpolation  of  the  values  Fk  =  {Fk}il0  at  the  grid  points  of  Xk  in 
(2.1a). 

In  the  previous  two  examples  we  have  shown  how  to  design  MR  schemes  for  discrete 
data  v ^  which  is  obtained  from  discretization  of  functions  by  point  values  (Example  1)  or  by 
cell-averages  (Example  2)  in  the  nested  sequence  of  uniform  grids  of  [0,1].  In  [H2]  and  [H3] 
we  have  presented  a  more  general  class  of  discrete  data  which  can  be  represented  by  a  scale 
decomposition.  This  class  is  characterized  by  the  following  notion  of  nested  discretization. 

Definition.  We  say  that  a  sequence  of  linear  operators  {Dk}^L0  is  a  nested  sequence  of 
discretization  if 

0) 

Vk:jr2^vk,  dim  Vk  =  Jk,  (2.22  a) 

(ii) 

Vkf  =  0  =*  T>k-if  =  0  (2.226) 

Here  T  is  a  space  of  mappings  and  Vk  is  a  linear  space  of  dimension  Jk. 

In  the  next  section  we  show  how  to  obtain  multiresolution  representation  of  any  discrete 
data  vL  =  VLf,  where  the  scale-decomposition  corresponds  to  the  levels  of  resolution  which 
are  introduced  in  (2.22).  This  is  a  very  general  framework  which  allows  for  discretizations 
corresponding  to  unstructured  grids  in  several  space  dimensions. 

3.  General  Multiresolution  Representation  Schemes 

In  this  section  we  consider  discrete  data  which  is  associated  with  a  nested  sequence  of 
discretization  {Dk}k=0  and  show  how  to  design  schemes  for  its  multiresolution  representation. 

First  we  show  that  a  nested  sequence  of  discretization  comes  equipped  with  a  decimation 
operator  Dkk~l  which  is  a  linear  mapping  from  Vk  =  Vk{F)  onto  Vk~x  =  Vk-\{T) 

Dk~'  :  Vk  ^  Vk~\  (3.1«) 

This  decimation  operator  is  defined  as  follows:  For  any  v  in  Vk  there  is  at  least  one  /  €  F 
such  that  Vkf  =  v ;  the  decimation  of  v  is  Vk-\f  €  Vk~x,  i.e., 

veVk,v=Vkf,  Dk~1v  =  Vk.1f.  (3.16) 
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It  follows  from  (2.22b)  that  Dk  1  is  well-defined  by  (3.1b),  i.e.  its  definition  is  independent 
of  the  particular  /.  To  see  that  let  us  take  fx  and  f2  in  T  such  that 

Vkfi  =  v  =  Vkf2 , 

then  by  (2.22b) 

0  =  X>fc/i  -  Vkf2  -  !>*(/,  -  /2)  =>  0  =  T>k-i(fi  -  /2)  =  XV-i/i  -  X>fc-i/2 

which  proves  our  claim. 

Given  vL  e  VL  we  can  evaluate  {v*}^1  by  repeated  decimation 

vk~l  =Dkk-xvh,  k  =  L,...,  1.  (3.2) 

Since  (3.1b)  implies  that 

Dk~l{T>kf)  =  Pfc-i/  for  any  /  €  T  (3.3) 

we  get  for  any  f  ef  for  which  vL  =  VLf,  that  u*  =  Vkf  for  all  k  in  (3.2).  We  would  like 
to  stress  the  point  that,  as  in  (2.2b)  and  (2.14b)  in  the  previous  examples,  this  decimation 
is  done  without  explicit  knowledge  of  /. 

Since  by  (2.22a)  Vk  =  Vk(T),  it  follows  that  Vk  has  a  right-inverse  (at  least  one)  which 
we  denote  by  Hk\ 

Kk  :Vk  T,  Vk7lk  =  Ik ,  (3.4) 

where  Ik  denotes  the  identity  operator  in  Vk .  Since  (' TZkvk )  G  T  is  an  approximation  to  any 
/  €  T  for  which  Vkf  —  vk ,  we  refer  to  TZk  as  a  reconstruction  of  Vk\  see  (2.3)-(2.4)  and 
(2.15)-(2.16)  in  the  previous  examples. 

Next  we  show  that  any  sequence  of  corresponding  reconstruction  operators  {7lk}k=0 
defines  a  MR  scheme  for  discrete  data  vL  in  VL.  Starting  from  in  (3.2)  we  can  get  an 
approximation  to  vk  by 

vk «  vk{nk.,vk~l). 

We  denote 

Pk~ i  =:  T)kTZk- 1,  Pk-\  '•  Vk  1  —y  Vk  (3.5a) 

and  refer  to  it  as  prediction  operator.  It  follows  immediately  from  taking  /  =  in 

(3.3)  and  using  (3.4)  that  P£_ j  is  a  right-inverse  of  the  decimation  Dkk~x 

Dk~'Pk-i  =  h-i-  (3.56) 

We  observe  that  the  prediction  error  ek 

ek  =  vk-  Pkk_xvk~x  =  (lk  -  Pk_1Dk~1)vk  (3.6a) 
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satisfies  the  relation 


Dkk~lek  =  Dkk~1vk  -  =  vk~k  -  vk~l  =  0 

and  therefore  it  is  in  the  null  space  of  the  decimation  operator 

ek€Af(Dkk~1)  =  {v\  veVk,Dkk~1v  =  0}. 


(3.66) 


(3.76) 


It  follows  from  (3.1a)  that 

dim-A/pJ"1)  =  Jk  ~  Jk-x  (3-7a) 

and  hence 

jV(^-1)  =  span{/.*}*7J‘-\  (3.76) 

where  is  any  basis  of  M(Dkk~l).  Therefore  the  prediction  error  ek,  which  is 

described  in  terms  of  Jk  components  in  Vk,  can  be  represented  by  its  (Jk  —  Jk- 1)  coordinates 
dk  in  (3.7b) 

ek  =  £  1  dknk  =:  Ekdk ,  dk  =:  Gkek.  (3.8a) 


(3.8a) 


Here  Gk  denotes  the  operator  which  assigns  to  ek  G  J\f  (D\  1 )  its  coordinates  dk  in  the  basis 
observe  that  EkGk  is  the  identity  operator  in  N(Dkk~x),  i.e. 

EkGkek  =  ek  for  any  ek  G  N(Dk~l).  (3.86) 

At  this  point  we  can  show  that  there  is  a  one-to-one  correspondence  between  vk  and 
{dfc,vfc-1}:  Given  vk  we  evaluate 

f  =  d*-V 

1  dk  =  Gk(h  -  PLiD^y  ; 
given  ufc_1  and  c/fc  we  recover  t>fc  by 

flX"1  +  Ekdk  =  PLllfr1v>  +  EkGk{Ik-Jfc.1Dir1V 
=  P£_1Dkk-1vk  +  (Ik  -  Pt-iDt'y 
=  vk. 

As  in  the  previous  examples  this  shows  that 


„L  1:1  (jL  Jl  .,0 


{dL,...,d1,u0}  =:  UAf, 
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where  the  direct  MR  transform  vm  =  M  ■  vL  is  given  by  the  algorithm 

DO  k  =  L, . . . ,  1 

u*"1  =  Dkk~1vk  (3.10) 

dk  =  Gk{Ik  -  PLxDk-')vk  =:  G?vk 
and  the  inverse  MR  transform  vL  =  M~l  •  vm  can  be  calculated  by 

DO  k  =  1 

(3.11) 

vk  =  p^y-1  +  Ekdk. 

We  remark  that  in  multigrid  terminology  Dkk~x  is  “restriction”  and  Pk-1  is  “prolongation.” 
In  signal  processing  Dkk~l  plays  the  role  of  “low-pass  filter”  while  Gk  ,  which  is  defined  in 
(3.10),  plays  the  role  of  “high-pass  filter.” 

Example  3.  Biorthogonal  Wavelets. 

In  this  example  we  derive  the  MR  schemes  which  correspond  to  the  bases  of  biorthogonal 
wavelets  in  [CDF].  These  MR  schemes  are  obtained  from  nested  discretization  of  functions  in 
^?oc(R)  by  taking  weighted- averages  on  a  nested  sequence  of  uniform  grids  of  R,  as  follows: 

1  f°°  ( x  — 

{Dkf)i  =  J  f(x)w  ( — - — LJ  dx,  — oo  <  i  <  oo,  (3.12a) 

where  w  €  T2(R)  is  a  weight-funcntion 


/*oo 

/  w(x)dx  =  1 

(3.12  b) 

and 

J  —  oo 

xk  =  ihk,  — oo  <  i  <  oo,  hk  =  2 ~khQ. 

(3.12 c) 

In  order  to  obtain  a  nested  sequence  of  discretization  we  want  to  choose  w(x)  so  that 

N 

{Dk-\f)i  =  e(Pkf)2i-t  (3.12d) 

e=o 

where  on ,  l  —  0, . . . ,  N  are  real  numbers.  We  observe  that  since 

f(x)  =  c=  constant  =>  ( T>kf)i  =  c  Vi,  k 
we  have  to  limit  the  choice  of  {«/■}  by 

N 

ae  =  1-  (3.13a) 

r=o 
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From  (3.12d)  and  (3.12a)  we  get  that  for  any  /  €  Ti2oc(R) 


this  shows  that  w(x)  has  to  satisfy  the  functional  equation 

N 

w(x )  =  2  a£w(  2x  —  £).  (3.136) 

e=o 

This  equation  has  been  investigated  in  [Da]  and  [CDF]  and  got  the  name  of  “dilation  equa¬ 
tion”  in  [S].  It  is  shown  there  that  subject  to  condition  (3.13a)  the  dilation  equation  (3.13b) 
has  a  distribution  solution  which  is  determined  up  to  a  multiplicative  constant  and  a  shift, 
and  that  w(x)  has  a  support  of  size  N.  Furthermore,  if 

<*2  e  —  E  a2i-i  (3.13c) 

e  i 

then  w(x)  is  also  square  integrable. 

We  make  the  choice  of  w(x)  unique  by  imposing  (3.12b)  and  fixing  its  support  in,  say, 
[— N,  0].  With  this  choice  of  a  weight-function,  the  sequence  of  discretization  in  (3.12)  is 
nested  and  its  decimation  operator  is  given  by 


(Dl  1v)i  =  ^2a2i-mvm.  (3.14a) 

t  m 

In  [Da]  and  [H3]  it  is  shown  that  0, 

(fij)i  =  (-l)’+1«2j-*-i,  -oo  <i  <oo  (3.146) 

is  a  basis  of  in  (3.6b). 

We  reconstruct  the  discretization  Djt  in  (3.12a)  as  follows:  We  take  a  sequence  {0e}  of 
compact  support  which  satisfies 


12e  fct  —  Yht  &2i-\  —  1 

i 

_  T)/  &£0£+2m  =  6m,o 


(3.15) 


and  define 


(ftfcu)(x)  = 


(3.16a) 


where  <p(x)  is  a  solution  of  the  dilation  equation 


¥>(*)  =  J2Pt<p{2x  -t) 

£ 


(3.166) 
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which  is  normalized  by 


J  <p(x)w(x)dx  =  1.  (3.16c) 

It  is  easy  to  see  that  the  corresponding  prediction  operator  P£_ j  =  'Dk'R-k- 1  is  given  by 

(Pk- (3.17) 

771 

Daubechies’  orthonormal  wavelets  are  obtained  by  imposing  the  additional  condition 

fit  =  2  at 

(see  [H2]  and  [H3]  for  more  details). 

We  remark  that  the  “fundamental  solution”  for  the  dilation  equation  (3.13b)  with 


<*t  =  f>t,0 

(3.18a) 

(x)  =  (5(x), 

(3.186) 

where  S(x )  is  the  Dirac  distribution  (see  [S] ) .  In  this  case  (3.12a)  becomes  the  point  value 
discretization  (2.2a)  of  Example  1.  However,  since  <$(x)  is  not  square  integrable,  point  value 
discretization  is  excluded  from  the  theory  of  wavelets.  For 

Q(  =  2^’°  (3.19a) 

we  get 

=  (  0  otherwise  (3'196) 

which  is  square  integrable,  and  the  discretization  (3.12a)  becomes  the  cell-average  discretiza¬ 
tion  (2.14a)  of  Example  2.  Observe,  however,  that  the  theory  of  wavelets  is  for  the  infinite 
domain  R,  while  our  formulation  is  suitable  for  both  the  finite  (Example  2)  and  the  infinite 
case  (Example  3).  Furthermore,  unlike  the  theory  of  wavelets  which  uses  translates  and 
dilates  of  a  single  function  for  both  discretization  and  reconstruction  and  consequently  is 
restricted  to  uniform  grids,  our  framework  of  nested  discretization  allows  for  general  ge¬ 
ometries.  In  [H3]  and  [AH1]  we  extend  the  multiresolution  representation  of  cell-averaged 
data  in  Example  2  to  unstructured  meshes  in  bounded  domains  of  Rm,m  >  1,  by  using 
agglomeration  of  cells  to  generate  a  nested  sequence  of  discretization. 

In  [ADH]  we  consider  the  case  where  w(x)  in  (3.12)  is  the  “hat-function”  and  show  how 
to  improve  data  compression  by  using  adaptive  prediction  techniques  near  discontinuities 
and  distributions. 
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4.  Multiresolution  Representation  of  a  2-Dimensional  Array 

In  this  section  we  consider  functions  / 

/  :  [0, 1]  x  [0, 1]  — R  (4.1a) 

which  are  discretized  on  the  tensor-product  grid 

=  (4.U) 

=  (Aij?  /„'  I'  f{x"x*)w  1^)  ”  {^rr)  dx'dx 2  (4-lc) 

where  {xf}  are  the  one- dimensional  gridpoints  in  (2.1)  and  w(x)  is  a  weight-function  as  in 
Example  3;  by  (3.18)  and  (3.19)  this  includes  point  value  and  cell-average  discretizations. 
Although  this  case  is  covered  by  the  general  framework  in  (3.9)-(3.11),  it  is  convenient  here 
to  represent  the  two-dimensional  array  in  (4.1)  as  the  Nk  x  Nk  matrix  Ak,  and  to  use  tensor- 
product  extension  of  the  corresponding  one-dimensional  operators  to  get  a  MR  scheme  for 
the  input  AL. 

Let  us  denote  the  matrix  representation  of  the  various  one- dimensional  operators  by 

Dkk~l  — ♦  (D)Nk_1xNk 

Pk- 1  - *  {P)NkxNk^ 

(4.2) 

Gk  =  Gk{Ik  -  PLM-1)  —  {CPfr-tx*  =  G(I  -  PD) 

Ek  - »  (E)NkxNk_1- 

These  matrix  representations  are  obtained  by  taking  vk  and  dk  in  (3.10)-(3.11)  to  be  column- 
vectors,  e.g. 

Nk 

v?-1  =  (Dt1  vk)i  '=:  ^2DijVk,  1  <  *  <  Nk _i  (4.3a) 

j=i 

where  by  (3.14a) 

Dij  =  a2i-j  ;  (4.36) 

for  simplicity  we  drop  the  index  k. 

Starting  with  AL  we  decimate  to  get 

Ak~1  =  DAkD*,  k  =  L,...,  1;  (4.4) 

here  (•)*  denotes  the  transpose.  Given  Afc_1  we  get  an  approximation  to  Ak  by 

Ak  ts  PAk~kP* 
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and  observe  that  the  prediciton  error  matrix  ek 


ek  =  Ak  -  PAk~lP*  (4.5a) 

satisfies 

DekD*  =  0.  (4.56) 

The  dimension  of  the  null  space  of  the  decimation  operators  here  is 

Jk  ~  Jk- 1  =  (2iVjt_i)2  -  (iVfc-i)2  =  3 (7V,_i)2 

and  we  store  the  scale-coefficients  dk  in  three  JV^-i  x  Nk- 1  matrices  A*,  A*,  A3. 

Using  the  matrix  identity 

I  =  PD  +  EGd,  Gd  =  G(I  -  PD),  (4.6) 

we  show  that  if  we  take  Ak  and  Ak-i  from  the  sequence  (4.4)  and  define 

Ak  =  GDAk(GDY ,  Ak  =  DAk{GDY ,  A*  =  GDAkD*  (4.7) 

then  Ak  can  be  recovered  from  Ak~ 1  and  the  above  by 

Ak  =  PAk~lP*  +  EA\E *  +  PAkE*  +  EAkP*.  (4.8) 


This  follows  immediately  from  the  identity 

Ak  =  (PD  +  EG°)Ak(PD  +  EGDy  =  P(DAkD*)P* 

+  E[GDAk(GDy)E*  +  P[DAk(GDy\E *  +  E(GDAkD*)P\ 


We  conclude  from  (4.7)-(4.8)  that  AL  has  a  multiresolution  representation  4^, 

A  (KL.--KL,-'40)’  <4-9> 

where  the  direct  MR  transform  is  given  by 
DO  k  =  L,...,l 

Ak~l  =  DAkD*  (4.10) 

A k  =  GDAk(GD)*,  A k  =  DAk(G°y,  Ak  =  G°AkD\ 

and  the  inverse  MR  transform  is 


{DO  k  = 

Ak  =  PAk~lP*  +  EA\E *  +  PA\E*  +  EAkP\ 


(4.11) 
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In  Figure  1,  which  is  taken  from  [HY],  we  show  the  result  of  applying  thresholding  with 
tolerance  e  =  10~7  to  Am  in  (4.9)  where 


t—-t  if  i  ^  j 
0  if  i  =  j 


(4.12) 


Here  we  take  Nl  —  512  and  use  point  value  discretization  which  is  reconstructed  by  a 
sixth-order  accurate  piecewise-polynomial  interpolation  with  a  centered  stencil,  near  the 
boundaries  we  use  one-sided  stencils.  The  rate  of  compression,  i.e.  the  ratio  of  {Nl)  to  the 
number  of  entries  in  A ^  which  are  above  the  tolerance  e  —  10  7,  is  8.57.  This  example  is 
taken  from  [BCR]  where  it  is  done  with  MR  schemes  which  use  Daubechies’  wavelets;  the 
corresponding  rate  of  compression  for  wavelets  with  six  vanishing  moments  is  reported  to  be 
7.33  . 

In  Figure  1  we  display  the  results  by  writing  Am  in  (4.9)  as  the  matrix 


and  marking  the  entries  that  are  larger  in  absolute  value  than  10  7  by  a  black  dot. 


5.  Multiresolution  Algorithm  for  Matrix-Vector  Multiplication. 

In  this  section  we  describe  an  algorithm  to  reduce  the  cost  of  evaluating  the  matrix-vector 
multiplication  c, 

c  =  Ab,  (5.1) 

where  A  is  a  Nl  x  Nl  matrix  which  can  be  thought  of  as  the  discretization  (4.1)  of  some 
piecewise-smooth  function,  and  b  is  any  vector  of  Nl  components;  note  that  we  do  not  make 
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any  smoothness  assumptions  on  6.  As  an  example  let  us  consider  the  case  where  (5.1)  is 
obtained  from  discretization  of  the  integral  transform 


:(x)  =  [  a(x,y)b(y)dy, 
Jo 


(5.2  a) 


here 

C{  ~  c(:ct-  )’  ^«i  ~  ’  o(®{  bj  fti  ),  (5.26) 

and  we  want  to  evaluate  (5.1)  to  a  specified  tolerance  which  is  taken  to  be  of  the  order  of 
the  discretization  error  in  (5.2). 

The  basic  idea  is  that  the  product  c  =  A6,  which  is  being  set  up  as  cL  =  A^b^,  has  a 
meaningful  analog 

ck  =  Akbk ,  k  =  L-  1,...,0  (5.3a) 


corresponding  to  the  fc-th  grid,  and  that 


ck  =  Pck  1  +  correction,  (5.36) 

where  P  is  the  matrix  representation  (4.2)  of  the  one-dimensional  prediction  operator,  and 
the  “correction”  comes  from  locations  in  [0, 1]  x  [0, 1]  where  a(x,y)  is  still  not  “sufficiently” 
resolved  on  the  ( k  —  l)-th  grid.  To  obtain  (5.3)  let  us  multiply  (4.8)  by  bk  to  get 

Akbk  =  PAk~1(P*bk)  +  EAk(E *bk)  +  PAk2(E*bk )  +  EAk3(P*bk ), 

from  which  we  see  that  if  we  define 


bk~l=P*bk ,  k  =  L,...,  1  (5.4a) 

and  denote 

sk~l  =  E*bk,  k  =  l,...,L  (5.46) 

then  we  can  rewrite  the  above  identity  as 

ck  =  Pck~l  +  E{  A^-1  +  A^-1)  +  P(Aksk~1).  (5.4c) 

Using  this  observation  we  get  the  following  algorithm  for  the  approximate  (up  to  a  specified 
error)  evaluation  of  the  matrix- vector  multiplication  (5.1): 

(A)  Preparation: 

(i)  Given  A,  set  AL  =  A  and  use  the  direct  MR  transform  (4.10)  to  obtain  the  MR 
representation  Am  (4.9). 
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(ii)  Apply  thresholding  to  Am  in  order  to  obtain  a  sparse  representation. 

(B)  Multiplication: 

Given  any  vector  6 


(i)  Set  bL  =  b 

DO  k  =  L,...,  1 
<  bk~l  =  P*bk 
k  sk~ 1  =  E*bk 

(ii)  Calculate  directly 

c°  =  A°b° 

(this  is  done  on  the  coarsest  grid). 

(iii) 

DO  for  k  =  1,... ,L 

< 

ck  =  p^M- 1  +  a  ksk-')  +  E(  Afs'5-1  +  A^-1) 


Set  c  =  cf. 


(5.5) 


(5.6a) 


(5.66) 


In  the  case  of  pointvalue  discretization  this  algorithm  turns  out  to  be  identical  to  the 
multilevel  matrix  multiplication  of  Brandt  and  Lubrecht  in  [BL].  The  algorithm  which  cor¬ 
responds  to  Daubechies’  orthonormal  wavelets  is  identical  to  the  “non-standard  form”  of 
Beylkin,  Coifman  and  Rokhlin  in  [BCR].  We  remark  that  (5.5)-(5.6)  is  a  slight  generaliza¬ 
tion  of  the  algorithm  which  was  presented  in  [HY]. 

Based  on  the  analysis  in  [BCR]  we  show  in  [HY]  that  if  the  kernel  a(x,  y )  in  (5.2a)  satisfies 

\dea{x,y)\  <  ^  (5-7) 

then  outside  a  diagonal  band  of  width  B,  the  entries  of  the  matrices  { Akri  }^1=j  are  not  larger 
than  0(B~r),  where  r  is  the  order  of  accuracy  of  the  reconstruction  technique.  Since  the 
matrices  P  and  E  are  banded,  we  get  that  the  complexity  of  the  algorithm  in  (5.5)-(5.6)  is 

0(Nl)  . 

In  [HY]  we  used  the  above  matrix-vector  multiplication  algorithm  to  apply  the  matrix 
in  (4.12)  to  a  vector  6  with  “randomly  generated”  components.  Using  single  precision  we 
obtained  for  the  case  in  Figure  1  a  relative  residual  error  ||A6  — c[|/]|6||  which  was  7.52  x  10-6 
in  the  l\  norm,  and  4.41  x  10-5  in  the  4©  norm. 
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6.  Multiresolution  Form  of  Numerical  Schemes. 

In  this  section  we  consider  the  numerical  algorithm 

vn+1  =  Avn  +  g  (6.1a) 

which  describes  either  an  iterative  procedure  or  a  numerical  scheme  for  the  solution  of  an 
initial  value  problem.  Taking  the  MR  transform  of  (6.1a)  we  get 

1  =  Mvn+1  =  (M  AM~1)(Mvn)  +  Mg  =:  AsvnM  +  gM]  (6.16) 


observe  that  /is  =  MAM -1,  the  multiresolution  form  of  the  matrix  (=  operator)  A,  is 
different  from  Am,  the  multiresolution  representation  (4.9),  where  A  is  treated  as  a  two- 
dimensional  array  and  not  as  an  operator. 

From  (3.10)  we  get  that  MvL  can  be  expressed  by 

'  ik  =  Gk(DUi---Dl-')-vL=-.BlvL,  k  =  l,...,L 
•  (6.2) 
v°  =  (D°1---D£~1)-vl  =:  Bglvl. 

From  (3.11)  we  get  that  M~1vm  can  be  expressed  by 


M-1vM  =  'Ec£dk  +  C£v0, 


(6.3a) 


where 


^  =  (^->•"04  k  =  l,...,L 


Using  (6.2)-(6.3)  to  epxress  (6.1b)  with  g  =  0,  we  get 

r  =  EL, (BlAC?)  ■  +  (BlACb)  ■  v°.”,  k  =  1, . . . ,  L 


,,0,71+1  _ 


Zf=t(B°LAC^)  ■  de’n  +  (B^ACq)  •  u0-" 


(6.36) 


(6.4a) 


which  shows  that 
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b°lac {■ 


Observe  that  the  block  B^ACf  ,  l<k,£<  L ,  is  of  size  Nk-i  x  Ne- 1. 

In  the  Appendix  we  show  that  Cf  =  (#£)*  where  is  given  by  (6.2)  for  the  dual  MR 
scheme;  hence 

(6.5)  BkLACi  =  BkLA(BlLy 

and  each  column  of  (BkLA)  is  dk.  the  scale  coefficients  (6.2)  of  the  A'-th  column  of  A,  while 
each  row  6f  A(BeL)*  is  (<?)*,  where  dl  is  the  column- vector  of  scale  coefficients  of  the  data  in 
the  £-th  row  of  A  which  is  obtained  by  the  dual  MR  scheme. 

We  remark  that  when  M  in  (6.1b)  corresponds  to  Daubechies’  orthonormal  wavelets,  As 
is  identical  to  the  standard  form  which  was  introduced  in  [BCR].  In  this  case  the  MR  scheme 
is  dual  to  itself,  i.e.  B £  =  B £  and  therefore  M-1  =  M*.  It  was  shown  in  [BCR]  that  if  A  is  a 
discretization  of  a  Calderon- Zygmund  operator,  then  data  compression  of  each  block  in  (6.5) 
results  in  a  “diagonal”  band,  the  width  of  which  is  independent  of  k  and  t.  Consequently 
applying  data  compression  to  As  results  in  a  finger-like  structure  of  non-zero  entries,  and 
their  total  number  is  0(Nl  log2  Nl)- 

In  [AC]  and  [ACD]  the  standard  form  of  [BCR]  was  extended  by  (6.1b)  to  point  value  and 
cell-average  discretizations.  It  is  shown  there  that,  inspite  of  the  lack  of  symmetry,  the  rate 
of  compression  compares  favorably  to  that  of  orthonormal  wavelets:  It  is  about  the  same  in 
the  periodic  case,  but  it  is  significantly  better  in  presence  of  boundaries. 

In  Figure  2,  which  is  taken  from  [AC],  we  show  the  nonzero  entries  of  As ,  the  multires- 
olution  form  of  the  matrix  in  (4.12).  As  in  Figure  1  we  use  point  value  discretization  which 
is  reconstructed  by  piecewise-polynomial  interpolation,  however  here  Nl  =  64. 
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The  present  work  extends  the  scope  of  application  of  these  algorithms  to  the  general  class 
of  MR  schemes  in  (3.9)-(3.11). 

7.  Multiresolution  Application  of  Operators. 

In  this  section  we  describe  multiresolution  algorithms  for  the  solution  of  integral  equations 
and  for  the  numerical  solution  of  initial-value  problems.  We  shall  say  that  a  matrix  is  T- 
Sparse  if  it  becomes  sparse  under  thresholding,  and  refer  to  the  number  of  non-zero  elements 
as  its  T-sparseness.  In  this  language  we  can  express  the  main  result  of  the  previous  section 
by  saying  that  the  multiresolution  form  As  of  Calderon- Zygmund  operators  is  T-sparse  and 
that  the  T-sparseness  of  As  is  0(Nl  log2  Nl). 

7. A.  Integral  Transforms  and  Equations. 

In  section  5  we  described  a  multilevel  algorithm  for  matrix- vector  multiplication  (5.1) 
and  showed  that  if  it  corresponds  to  a  discretization  of  an  integral  transform  of  the  form 
(5.2a)  with  a  kernel  which  satisfies  (5.7),  then  it  can  be  performed  in  0(NL)  operations.  We 
can  also  evaluate  c  =  Ab  from  its  multiresolution  form  by 

'  bM  =  Mb, 

i  cm  =  AsI>m ,  (7.1) 

.  c  = 

Since  the  T-sparseness  of  As  is  0(NL  log2  NL)  this  is  also  the  cost  of  the  product  AsbM\  in 
addition  we  also  have  to  apply  the  direct  MR  transform  to  b  and  its  inverse  to  cjv/.  Hence, 
unless  we  are  in  a  special  situation  in  which  cm  and/or  b\i  are  also  T-sparse,  it  is  more 
efficient  to  calculate  integral  transforms  by  the  multilevel  algorithm  (5.5)-(5.6). 

A  situation  of  this  type  occurs  in  a  matrix-matrix  multiplication  C  =  AB  where  both 
As  and  Bs  are  T-sparse 

Cs  =  MCM-1  =  (MAM~1)(MBM~1)  =  ASBS.  (7.2a) 

Of  particular  interest  is  the  case  where  the  T-sparseness  of 

M(A)nM~l  =  (As)n  (7.2  b) 

is  uniform  in  n.  In  this  case  (As)n  for  n  —  2m  can  be  computed  in  m  steps  of  squaring  and 
thresholding 

(As)2"  =  [(As)2*-1]2,  k  =  1, . . .  ,m,  (7.3) 

where  each  product  above  is  between  matrices  with  T-sparseness  of  0(NL  log2  NL). 
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The  numerical  algorithm  (6.1a)  can  be  written  as 

u"=(A)v+£W~y  (7-4«) 

j= 0 

or  in  its  multiresolution  form  (6.1b)  as 

vnM  =  (As)nv°M  +  Yl(Asy-1gM.  (7.46) 

j=o 

Following  [BCR],  [EOZ]  and  [ACD]  we  use  the  following  algorithm  for  the  fast  evaluation  of 
(7.4a)  for  n  =  2m: 


(i)  Set 


B  =  MAM-1 


(ii) 


Calculate 


C  =  I 

DO  m  times 
C  =  tr(C  +  BC-,e) 
B  =  tr(BB-,e ) 


vn  =  M~1(BMv°  +  CMg); 
here  tr(A ;  e)  denotes  the  truncation  operation 


(7.5) 


.  '  Aij  if  \AiJ  >e 

[tr(A]  e)]y  =  <  .  (7.6) 

.  0  if  I  A{j  |  <  e 

Fredholm  integral  equations  of  the  second  kind  are  usually  solved  by  iterative  procedures 
of  the  form  (6.1a).  In  this  case  the  gain  in  efficiency  which  is  offered  by  algorithm  (7.5)  is 
based  both  on  the  compression  of  the  operator  A  by 

A  — >  tr(A3-,e)  (7.7) 

and  using  the  fact  that  the  T-sparseness  of  (As)n  remains  constant  as  n  increases. 
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7.B.  Initial  Value  Problems. 

Consider  the  evolution  equation 

dtu  +  C(x,dx)u  =  f(x)  <>  0  (7  . 

u(x,0)  =  uo(ar)  1  > 

with  boundary  conditions,  where  £  is  a  differential  operator.  An  explicit  discretization 
typically  takes  the  form  (6.1a) 

vn+ 1  =Avn  +  g  (7.9a) 

where 

Vj  «  u(xj,tn ),  g3  fa  Atf(xj),  (7.9 6) 

tn  =  nAt  and  {x.,}  is  a  grid  in  Cl. 

It  is  shown  in  [EOZ]  that  using  the  multiresolution  algorithm  (7.5)  to  calculate  large¬ 
time  solutions  of  one-dimensional  hyperbolic  problems  to  a  fixed  predetermined  accuracy 
can  be  reduced  from  the  standard  0((Nl)2)  to  0(AT(log2  AT)3).  For  parabolic  equations, 
a  standard  explicit  calculation  with  complexity  0((Nl)3)  can  be  likewise  reduced  by  (7.5) 
to  0(AT(log2  AT)3)-  The  multiresolution  algorithm  (7.5)  in  [EOZ]  is  based  on  Daubechies’ 
orthonormal  wavelets.  In  [ACD]  this  algorithm  is  extended  to  point  value  and  cell-average 
discretizations. 

As  an  example  let  us  consider  the  simplest  hyperbolic  problem 


Uf  -j-  U%  0 


(7.10a) 


and  its  solution  by  the  Lax-Wendroff  scheme 


u,n+1  =  v? 


-  »?-,)  +  -  2»r  +  »?-,)  =:  (^n)i 


(7.106) 


where  A  =  At/h^. 

The  matrix  A  is  a  tridiagonal  matrix  and  thus  has  3 Ni  non-zero  entries.  On  the 
other  hand  As,  the  multiresolution  form  of  the  scheme,  has  a  finger-like  structure  with 
0( Nl  log2  Nl)  non-zero  entries.  In  Figure  3,  which  is  taken  from  [AC],  we  show  the  mul¬ 
tiresolution  form  of  the  Lax-Webdroff  scheme  for  Nl  =  64.  This  shows  that  unlike  the 
application  to  iterative  solution  of  integral  equations  where  (7.7)  results  in  a  “compressed” 
representation  of  the  operator,  the  gain  in  efficiency  in  large-time  computation  of  hyperbolic 
problems  is  only  due  to  the  uniform  T-sparseness  of  powers  of  As,  i.e.  while  (A)n  fills  up 
linearly  in  n,  the  T-sparseness  of  (As)n  remains  0(Nl  log2  AT)  for  all  n. 

In  the  following  we  describe  another  application  of  the  multiresolution  form  (6.1b)  which 
uses  data-compression  of  the  numerical  solution  vn  and  the  finite-speed  of  propagation  in 
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hyperbolic  problems  in  order  to  produce  a  multiresolution  analog  to  adaptive  grids.  Let 

re(0  =  {(j,*)l  KV)I>£*>  (7-n«) 

denote  the  domain  in  the  ( j  -  k )  plane  which  contains  all  the  significant  scale-coefficients 
of  vn,  and  let  f”  denote  the  domain  which  is  obtained  by  enlarging  re(un)  by  adding  side- 
neighbors  of  the  cells  in  re(u”)  and  allowing  for  a  growth  of  one  scale  per  time-step  where 
needed.  Due  to  the  finite  speed  of  propagation  in  hyperbolic  problems 

rE(un+1)  C  f?.  (7.116) 

Therefore  we  can  set  the  components  of  which  are  not  listed  in  T”  to  zero,  and  evaluate 
the  product  row(As)  times  only  for  those  rows  which  are  listed  in  f  ”.  The  computational 
work  can  be  further  reduced  by  taking  into  account  the  T-sparseness  of  v 

This  technique  can  be  extended  to  nonlinear  problems.  In  [LT],  [MR]  and  [BMP]  it  is 
shown  how  to  derive  a  multiresolution  scheme  for  the  numerical  solution  of  Burgers’  equation 

Uf  T  wHx  —  vUxxi  v  0  (7.12) 

in  which  the  time-evolution  is  restricted  to  the  significant  scale-coefficients  of  the  numerical 
solution.  This  numerical  scheme  is  obtained  by  a  Galerkin  approach  in  which  the  PDE  is 
projected  on  a  basis  of  wavelets.  We  remark  that  this  Galerkin-type  scheme  is  not  suitable 
for  the  “inviscid”  Burgers’  equation  (u  =  0  in  (7.12))  in  the  sense  that  it  generates  spurious 
oscillations  at  shocks,  and  may  even  become  unstable  in  some  cases  -  thus  some  form  of 
artificial  viscosity  is  needed. 

In  [H4]  and  [H5]  we  consider  a  hyperbolic  system  of  conservation  laws 

««  +  /(«)*  =  0  (7.13a) 

and  its  numerical  solution  by  any  scheme  in  conservation  form 

vj+1  =  v;  -  \(fi  -  (7.13ft) 

where  ' 

fj  =  f(vj~K+ 1,  •  •  • ,  v^+K)  for  some  K  >  1,  (7.13c) 

and  /  is  the  numerical  flux  function.  Observe  that  the  computational  task  here  is  the  evalu¬ 
ation  of  the  numerical  flux  function  (7.13c)  at  all  the  gridpoints.  Using  the  multiresolution 
form  of  the  numerical  scheme  (7.13b)  with  respect  to  cell-average  discretization  we  derive  an 
algorithm  for  the  time-evolution  of  the  scale-coefficients,  and  show  that  data  compression  of 
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the  numerical  solution  can  be  translated  into  reduction  of  the  number  of  flux  calculations  in 
(7.13c). 

In  Figure  4  we  show  the  results  of  [H5]  for  the  multiresolution  form  of  the  Lax-Wendroff 
scheme  which  is  applied  to  the  periodic  initial  value  problem  for  the  “inviscid”  Burgers’ 
equation. 

ut  +  uux  =  0  ,  —  1  <  x  <  1  ,t  >  0 


u(x,0)  =  2  +  sin(7rx) 


(7.14) 


with  Jl  =  256,  CFL  =  0.8,  and  tolerance  e  =  10-3.  Figure  4  consists  of  3  snapshots 
corresponding  to  n  =  25, 150, 400  time-steps.  In  the  upper  part  of  each  snapshot  we  compare 
the  solution  of  the  MR  scheme  (circles)  to  the  solution  of  Lax-Wendroff  scheme  on  the  finest 
grid  (continuous  line),  which  is  computed  independently.  In  the  lower  part  of  each  snapshot 
we  display  rs(un)  (circles)  and  its  corresponding  estimate  f"  (dots).  This  is  done  by  drawing 
the  symbol  at  (a^J-i?  k)  in  the  x  —  k  plane  for  each  (j,  k )  in  the  set;  note  that  due  to  a  different 
notation  in  [H5]  k  =  0  is  the  finest  grid,  and  k  =  L  =  5  is  the  coarsest.  In  the  Table  we  list 
the  efficiency  (i.e.  the  ratio  between  the  fine  grid  calculation  of  256  fluxes  over  the  number 
that  we  actually  had  to  compute)  and  the  difference  Em,m  =  1,2,  oo  in  the  corresponding 
norm  between  the  solution  of  the  MR  scheme  and  the  independent  finest-grid  calculation. 

In  [BH1]  we  extend  this  technique  to  the  numerical  solution  of 


ut  +  divf(u)  =  0  (7.15) 

on  Cartesian  grids,  and  in  [AH2]  we  generalize  it  further  to  unstructured  meshes  where  the 
coarser  levels  of  resolution  are  obtained  by  agglomeration  of  cells. 


8.  Conclusions. 

In  this  paper  we  reviewed  recent  developments  in  techniques  to  represent  data  in  terms 
of  its  local  scale  components.  These  techniques  enabled  us  to  obtain  data  compression  by 
eliminating  scale- coefficients  which  are  sufficiently  small.  This  capability  for  data  compres¬ 
sion  can  be  used  to  reduce  the  cost  of  many  numerical  solution  algorithms  by  either  applying 
it  to  the  numerical  solution  operator  in  order  to  get  an  approximate  sparse  representation, 
or  by  applying  it  to  the  numerical  solution  itself  in  order  to  reduce  the  number  of  quantities 
that  need  to  be  computed. 
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Appendix:  The  Dual  MR  Scheme. 

In  this  appendix  we  describe  the  MR  scheme  which  is  dual  to  the  one  in  (3.10)-(3.11). 
In  order  to  better  see  the  duality  we  rewrite  (3.10)-(3.11)  as  follows:  we  express  the  direct 
MR  transform 

i>M  —  M  •  v ^  (A. la) 

by 

DO  k  = 

i  vk~l  =  Dk~tvk  (A.  16) 

dk  =  G%vk 

where 

Gf  =  Gk(h  -  Pl.Dt-').  (/1.1c) 

We  observe  that 

Ekdk  €  M(Dk~l)  =>  (Ik  -  Pk_1Dk~1)Ekdk  —  Ekdk 
and  therefore  we  can  rewrite  the  inverse  MR  transform 

vL  =  M~1vm  (A. 2a) 


f  DO  k  =  1, . . . , 
|  vk  =  Pk_xvk~l  - 


,L 

+  Ekdk 


(A.26) 


where 


Ek  =  (A  -  pk-iDk~')Ek-  (A 2c) 

To  simplify  our  presentation  we  shall  use  the  matrix  representation  of  the  various  operators 
and  define 

Dt' =:  (Pl  ,)*,  Plk=-(Dl~')\  Gf=:(£fr,  £f=:(Gf)*.  (A3) 

Observe  that  Dk _1  is  a  Jk- 1  x  Jk  matrix,  Pk_x  is  a  Jk  x  Jfc_i  matrix,  Gk  is  a  ( Jk  -  Jk- 1)  x  Jk 
matrix,  and  Ek  is  a  Jk  x  ( Jk  —  Jk-\)-  With  these  definitions  we  obtain  the  dual  MR  scheme 
from  (A.1)-(A.2)  by  applying  (•)  to  all  the  operators,  i.e.  the  direct  MR  transform  of  the 
dual  scheme 

va=MvL  (/1.4c) 


is  given  by 


DO  k  =  L, . . . ,  1 


vk  1  =  Dk  1vk 


(A.46) 


dk  =  Gf  vk 
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where  Gjf  is  defined  in  (A. 3).  The  inverse  MR  transform  of  the  dual  scheme 

vL  =  M-1  •  Vfo  (A. 5a) 

is  given  by 

'  DO  k  =  l,...,L 

<  (A. 56) 

.  vk  =  Pty-1  +  Efdk 

where  E£  is  defined  in  (A. 3).  Observe  that  the  dual  of  the  dual  is  the  original  scheme. 

It  follows  from  the  above  definitions  that  for  1  <  k  <  L 


ct  =  (*£.,  -n‘+,W  =  KEfn^+'r  •••(/>/■_,)•]•  =  ■■••Dir1)]-  =  (Sir 

and 


Go  =  Pt-1  -P^  [(n1)* •  •  •  (Pl-,)'Y  =  (Dl-- £)£-■)]-  =  (Sir  i 


thus 


AT1  =  {My  . 


In  [H3]  we  also  show  that  M  is  associated  with  discretization  T>k  and  reconstruction  7?* 
such  that  (TZkVk)  :  T  — >•  T  is  the  adjoint  of  (' RkVk )  . 
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Figure  1.  Multiresolution  representation  of  the  2-D  array  . 


Figure  2.  Multiresolution  form  of  the  operator  (matrix). 
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Figure  4.  Multiresolution  application  of  the  Lax-Wendroff  scheme. 
(Inviscid  Burgers’  equation). 
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